Spatial trajectory

11/18/24

Load data and libraries

##################
# LOAD LIBRARIES #
##################
library(tidyverse)
library(Seurat)
library(SeuratObject)
library(tidyseurat)
library(cowplot)
library(patchwork)
library(openxlsx)
library(slingshot)

source("../bin/spatial_visualization.R")
source("../bin/plotting_functions.R")

#########
# PATHS #
#########
input_dir <- "../results/04_deconvolution_st_data/"
result_dir <- "../results/08_spatial_dist/"
if( isFALSE(dir.exists(result_dir)) ) { dir.create(result_dir,recursive = TRUE) }

#############
# LOAD DATA #
#############
# DEGs_table <- read_csv(paste0(input_dir,"subset_100/DGEs_condition_wilcox_epi_100.csv")) %>% filter(groups == "DMPA")
DEGs_table_cond <- read_csv(paste0("../results/06_DGE_condition_st_data/","DGEs_condition_wilcox.csv"))
DEGs_table_clus <- read_csv(paste0("../results/05_DGE_clusters_st_data/","DGEs_clusters_wilcox.csv"))

DATA <- readRDS(paste0(input_dir,"seuratObj_deconvolution_scdc.RDS"))

#################
# COLOUR PALLET #
#################
col_epi <- c("#E41A1C","#FF7F00","#C77CFF","#984EA3")
col_submuc <- c("#CD9600","#7CAE00","#e0e067","#00A9FF","#377EB8","#984EA3",NA, NA,"#00BFC4", NA,"#FFFF33")
col_trajectory <- c("grey70", "orange3", "firebrick", "purple4")

col_feat <- c("#EFEDF5", "#DADAEB", "#BCBDDC", "#9E9AC8", "#807DBA", "#6A51A3", "#54278F", "#3F007D") # Purples
ord <-  c("Superficial","Upper IM","Lower IM","Basal", "8","3","4","0","2","1")
ID <- c("P107", "P108", "P114", "P097","P118", "P105", "P080", "P031")

Calculate spatial distance based on gene expression

DATA <- RunUMAP(DATA,
                n.components = 3L,
                dims = 1:50, 
                n.neighbors = 10,
                min.dist = .1,
                spread = 1,
                repulsion.strength = 1,
                negative.sample.rate = 10,
                n.epochs = 100,
                reduction = "harmony",
                reduction.name = "umap3Dharmony"
                )
pseudodist.fun <- function(DATA, start_clus, end_clus, Lineage="Lineage1"){
  # Define lineage ends
  ENDS <- c(end_clus)
  
  set.seed(1)
  lineages <- as.SlingshotDataSet(getLineages(
    data           = DATA@reductions$umap3Dharmony@cell.embeddings,
    clusterLabels  = DATA$Clusters,
    dist.method    = "mnn",       # It can be: "simple", "scaled.full", "scaled.diag", "slingshot" or "mnn"
    end.clus       = end_clus,    # You can also define the ENDS!
    start.clus     = start_clus)) # define where to START the trajectories
  
  
  # IF NEEDED, ONE CAN ALSO MANULALLY EDIT THE LINEAGES, FOR EXAMPLE:
  # sel <- sapply( lineages@lineages, function(x){rev(x)[1]} ) %in% ENDS
  # lineages@lineages <- lineages@lineages[ sel ]
  # names(lineages@lineages) <- paste0("Lineage",1:length(lineages@lineages))
  # lineages
  
  
  # Change the reduction to our "fixed" UMAP2d (FOR VISUALISATION ONLY)
  lineages@reducedDim <- DATA@reductions$umap_harmony@cell.embeddings
  
  # Calculate cluster centroids (for plotting the labels later)
  mm <- Matrix::sparse.model.matrix(~0 + factor(DATA$Clusters))
  colnames(mm) <- levels(factor(DATA$Clusters))
  centroids2d <- as.matrix(t(t(DATA@reductions$umap_harmony@cell.embeddings) %*% mm)/Matrix::colSums(mm))
  
  # plot
  # pal <- c((scales::hue_pal())(8), RColorBrewer::brewer.pal(9, "Set1"), RColorBrewer::brewer.pal(8,"Set2"))
  # set.seed(1)
  # pal <- rep(sample(pal, length(pal)), 200)
  # plot(DATA@reductions$umap_harmony@cell.embeddings, col = pal[DATA$Clusters],  cex=.5,pch = 16)
  # lines(lineages, lwd = 1, col = 'black', cex=2 )
  # text(centroids2d, labels = rownames(centroids2d),cex=0.8,font=2,col = "white")
  
  # Define curves
  curves <- as.SlingshotDataSet(getCurves(data = lineages, thresh = 0.1, stretch = 0.1,
      allow.breaks = F, approx_points = 300))

  print(curves)
  assign(paste0("line_", DATA$sp_annot[1]), curves@curves[[Lineage]][["s"]], envir = .GlobalEnv)
  
  pseudotime <- slingPseudotime(curves, na = FALSE)
  cellWeights <- slingCurveWeights(curves)
  
  x <- rowMeans(pseudotime)
  x <- x/max(x)
  
  name = paste0("sp_dist_", DATA$sp_annot[1])
  pseudotime <- pseudotime %>% 
    as_tibble(rownames = "barcode") %>% 
    dplyr::rename(!!name:=all_of(Lineage)) %>%
    select(-contains("Lineage")) %>%
    mutate(col = x) %>%
    arrange(desc(col))
    
  return(pseudotime)
}
epi_clus <- "^5$|^6$|^7|^9"
sub_clus <- "^5$|^8$|^4|^3|^2|^1$|^0$"

DATA_epi <- filter(DATA, grepl("epi", DATA$sp_annot) & grepl(epi_clus, DATA$Clusters) )
DATA_SubMuc <- filter(DATA, grepl("SubMuc", DATA$sp_annot) | DATA$Clusters == "5")
DATA_SubMuc <- filter(DATA, grepl(sub_clus, DATA$Clusters)| grepl("SubMuc", DATA$sp_annot) & DATA$sp_annot == "10" )

epi <- pseudodist.fun(DATA_epi, start_clus="6", end_clus="5")
class: SlingshotDataSet 

 Samples Dimensions
    1854          2

lineages: 1 
Lineage1: 6  9  7  5  

curves: 1 
Curve1: Length: 18.042  Samples: 1854
submuc <- pseudodist.fun(DATA_SubMuc, start_clus="5", end_clus="4")
class: SlingshotDataSet 

 Samples Dimensions
    5110          2

lineages: 1 
Lineage1: 5  8  3  0  1  2  4  

curves: 1 
Curve1: Length: 25.418  Samples: 5110
spatial_dist <- DATA %>%
  select(barcode=".cell" ) %>%
  left_join(., epi[,-3], by="barcode") %>%
  left_join(., submuc[,-3], by="barcode")
line <- line_epi
centroids2d <- DATA_epi %>%
    group_by(Clusters) %>%
    select(Clusters, matches("umapharmony_\\d", perl=T)) %>% 
    summarize_all(mean) %>% 
    column_to_rownames(var = "Clusters")

# dev.new(width=6, height=3, noRStudioGD = TRUE)
# Plot curves on UMAP coloured by clusters
par(mfrow = c(1,2), mar = c(3,3,0.1,0.1))
plot(DATA_epi@reductions$umap_harmony@cell.embeddings, 
     col = col_epi[DATA_epi$layers], pch = 16, cex = .5,  xlab = "", ylab = "",)
lines(line, lwd = 2, col = "black")
text(centroids2d, labels = rownames(centroids2d), cex = 1, font = 1)
title(ylab="UMAP 2", line=1.7, cex.lab=.9)
title(xlab="UMAP 1", line=1.7, cex.lab=.9)

plot(DATA_epi@reductions$umap_harmony@cell.embeddings[epi$barcode, ], 
     pch = 16,  xlab = "", ylab = "", # main = paste0("pseudotime"), axes = F 
     cex = 0.4, col = colorRampPalette(col_trajectory)(99)[epi$col * 98 + 1])
points(centroids2d, cex = 2.5, pch = 16, col = "#FFFFFF99")
title(ylab="UMAP 2", line=1.7, cex.lab=.9)
title(xlab="UMAP 1", line=1.7, cex.lab=.9)
text(centroids2d, labels = rownames(centroids2d), cex = 1, font = 1)

line <- line_SubMuc
centroids2d <- DATA_SubMuc %>%
    group_by(Clusters) %>%
    select(Clusters, matches("umapharmony_\\d", perl=T)) %>% 
    summarize_all(mean) %>% 
    column_to_rownames(var = "Clusters")

# dev.new(width=6, height=3, noRStudioGD = TRUE)
# Plot curves on UMAP coloured by clusters
par(mfrow = c(1,2), mar = c(3,3,0.1,0.1))
plot(DATA_SubMuc@reductions$umap_harmony@cell.embeddings, 
     col = col_submuc[DATA_SubMuc$Clusters], pch = 16, cex = .5,  xlab = "", ylab = "",)
lines(line, lwd = 2, col = "black")
text(centroids2d, labels = rownames(centroids2d), cex = 1, font = 1)
title(ylab="UMAP 2", line=1.7, cex.lab=.9)
title(xlab="UMAP 1", line=1.7, cex.lab=.9)

plot(DATA_SubMuc@reductions$umap_harmony@cell.embeddings[submuc$barcode, ], 
     pch = 16,  xlab = "", ylab = "", # main = paste0("pseudotime"), axes = F 
     cex = 0.4, col = colorRampPalette(col_trajectory)(99)[submuc$col * 98 + 1])
points(centroids2d, cex = 2.5, pch = 16, col = "#FFFFFF99")
title(ylab="UMAP 2", line=1.7, cex.lab=.9)
title(xlab="UMAP 1", line=1.7, cex.lab=.9)
text(centroids2d, labels = rownames(centroids2d), cex = 1, font = 1)

Save seurat object

DATA <- DATA %>%
  select(-contains("sp_dist")) %>%
  left_join(., spatial_dist, by=c(".cell"="barcode")) %>%
  select(c(1:7), contains("sp_dist"), everything())

saveRDS(DATA, paste0(result_dir,"seuratObj_spatial_dist.RDS"))
# DATA <- readRDS(paste0(result_dir,"seuratObj_spatial_dist.RDS"))

Plot spatial distance on tissue

####################
# SPATIAL DISTANCE #
####################
# dev.new(width=9, height=3.5, noRStudioGD = TRUE)
DATA %>%
  plot_st_feat.fun(.,
                   geneid = "sp_dist_epi",
                   zoom = "zoom",
                   col = col_trajectory,
                   alpha = 1,
                   ncol = 4, 
                   scale = FALSE,
                   annot_line = .1,
                   img_alpha = 0,
                   point_size = 1)

DATA %>%
  mutate(sp_dist_SubMuc = case_match(
    Clusters, "5" ~ NA, .default = sp_dist_SubMuc)) %>%
  plot_st_feat.fun(.,
                   geneid = "sp_dist_SubMuc",
                   zoom = "zoom",
                   col = col_trajectory,
                   alpha = 1,
                   ncol = 4, 
                   scale = FALSE,
                   annot_line = .1,
                   img_alpha = 0,
                   point_size = 1)

###########################
# CONDITION DEGs TO PLOT #
###########################
(sig_nest <- DEGs_table_cond %>%
  mutate(Direction = ifelse(avg_log2FC > 0, "UP", "DOWN")) %>%
  filter(p_val_adj < 0.05) %>%
  nest(., .by = c("layers", "Direction")) %>% 
  mutate(n = map_dbl(data, nrow)) %>%
  arrange(match(layers, ord), Direction))
# A tibble: 20 × 4
   layers      Direction data                   n
   <chr>       <chr>     <list>             <dbl>
 1 Superficial DOWN      <tibble [12 × 13]>    12
 2 Superficial UP        <tibble [7 × 13]>      7
 3 Upper IM    DOWN      <tibble [6 × 13]>      6
 4 Upper IM    UP        <tibble [11 × 13]>    11
 5 Lower IM    DOWN      <tibble [8 × 13]>      8
 6 Lower IM    UP        <tibble [16 × 13]>    16
 7 Basal       DOWN      <tibble [5 × 13]>      5
 8 Basal       UP        <tibble [11 × 13]>    11
 9 8           DOWN      <tibble [28 × 13]>    28
10 8           UP        <tibble [26 × 13]>    26
11 3           DOWN      <tibble [12 × 13]>    12
12 3           UP        <tibble [15 × 13]>    15
13 4           DOWN      <tibble [10 × 13]>    10
14 4           UP        <tibble [7 × 13]>      7
15 0           DOWN      <tibble [8 × 13]>      8
16 0           UP        <tibble [12 × 13]>    12
17 2           DOWN      <tibble [12 × 13]>    12
18 2           UP        <tibble [57 × 13]>    57
19 1           DOWN      <tibble [3 × 13]>      3
20 1           UP        <tibble [13 × 13]>    13
# all unique significant epi genes:
cond_epi_DEGs <- sig_nest %>% filter(!(grepl("\\d", .$layers))) %>% unnest(cols = c(data)) %>% pull(., "gene") %>% unique()
# all unique significant SubMuc genes:
cond_SubMuc_DEGs <- sig_nest %>% filter((grepl("\\d", .$layers))) %>% unnest(cols = c(data)) %>% pull(., "gene") %>% unique()

feat_cond <- c(cond_epi_DEGs, cond_SubMuc_DEGs)

###########################
# CLUSTER DEGs TO PLOT #
###########################
(sig_nest <- DEGs_table_clus %>%
  mutate(Direction = ifelse(avg_log2FC > 0, "UP", "DOWN")) %>%
  filter(p_val_adj < 0.05) %>%
  group_by(Clusters, Direction) %>%
  top_n(15, abs(avg_log2FC)) %>% 
  nest() %>% mutate(n = map_dbl(data, nrow)) %>%
  arrange(match(Clusters, ord), Direction) %>%
  ungroup())
# A tibble: 22 × 4
   Clusters    Direction data                   n
   <chr>       <chr>     <list>             <dbl>
 1 Superficial DOWN      <tibble [15 × 10]>    15
 2 Superficial UP        <tibble [15 × 10]>    15
 3 Upper IM    DOWN      <tibble [15 × 10]>    15
 4 Upper IM    UP        <tibble [15 × 10]>    15
 5 Lower IM    DOWN      <tibble [15 × 10]>    15
 6 Lower IM    UP        <tibble [15 × 10]>    15
 7 Basal       DOWN      <tibble [15 × 10]>    15
 8 Basal       UP        <tibble [15 × 10]>    15
 9 8           DOWN      <tibble [15 × 10]>    15
10 8           UP        <tibble [15 × 10]>    15
# ℹ 12 more rows
# all unique significant epi genes:
clus_epi_DEGs <- sig_nest %>% filter(!(grepl("\\d", .$Clusters))) %>% unnest(cols = c(data)) %>% pull(., "gene") %>% unique()
# all unique significant SubMuc genes:
clus_SubMuc_DEGs <- sig_nest %>% filter((grepl("\\d", .$Clusters))) %>% unnest(cols = c(data)) %>% pull(., "gene") %>% unique()

feat_clus <- c(clus_epi_DEGs, clus_SubMuc_DEGs)

############
# OVERLAP #
############
int <- intersect(feat_cond, feat_clus)
int
 [1] "SPRR2E"    "SPRR2A"    "SPRR2D"    "KRT78"     "KRTDAP"    "SPRR1B"   
 [7] "S100A7"    "IGHA1"     "JCHAIN"    "CSTA"      "IL1RN"     "IGHG1"    
[13] "IGHA2"     "KRT6A"     "PKP1"      "CRCT1"     "KRT14"     "LUM"      
[19] "KRT5"      "MT1X"      "CTSK"      "FLG"       "SPINK5"    "NEAT1"    
[25] "SERPINB3"  "KRT4"      "AQP3"      "TMPRSS11B" "REV3L"     "LGALS7"   
[31] "IGFBP5"    "MEG3"      "S100A9"    "SLC18A2"   "RPS27A"    "SFRP4"    
[37] "TMSB4X"    "CCN5"      "ATP5F1E"  
# select first gene only to run the script through fast:
cond_epi_DEGs <- cond_epi_DEGs[1]
cond_SubMuc_DEGs <- cond_SubMuc_DEGs[1]
clus_epi_DEGs <- clus_epi_DEGs[1]
clus_SubMuc_DEGs <- clus_SubMuc_DEGs[1]

Plot condition diff gene expression as dotplot and on tissue

##############################
# COND DOTPLOT PER EPI LAYER #
##############################
layer_dotplot.fun <- function(DATA, feat, spatial_dist, facet = TRUE, x_max=NULL, morf="epi", clus="^5$|^6$|^7|^9"){
  DAT <- DATA %>%
    filter(., grepl(morf, .$sp_annot)) %>%
    filter(., grepl(clus, .$Clusters)) %>%
    mutate(., FetchData(., vars = c(feat)) ) %>%
    select(orig.ident, groups, layers, all_of(c(feat)), {{spatial_dist}})
  
  if(morf=="epi"){probs <- c(0.009, 0.9)}else{probs <- c(0.1, 0.9)}
  
    rects <- DAT %>%
    group_by(layers) %>%
    summarise(., ystart=min({{spatial_dist}}, na.rm=T), yend=max({{spatial_dist}}, na.rm=T),
              Q1=quantile({{spatial_dist}}, probs = probs[1], na.rm=T),
              Q3=quantile({{spatial_dist}}, probs = probs[2], na.rm=T)) %>%
    filter(!(is.infinite(.$ystart))) %>%
    mutate(Q1 = ifelse(.$Q1 == min(.$Q1), 0,.$Q1)) %>%
    mutate(Q3 = ifelse(.$Q3 == max(.$Q3), max(.$yend),.$Q3)) %>%
    mutate(Q3 = ifelse(.$layers == "Lower IM", .$Q3+.3,.$Q3)) %>%
    arrange(ystart) %>% ungroup()

  if(facet == TRUE){facets <- facet_wrap(~groups, ncol = 2) }else{facets <- NULL}
  
  dot <- ggplot() +
    #ggtitle(feature) +
    geom_rect(data = rects, alpha = 0.1, show.legend=FALSE,
              aes(xmin = -Inf, xmax = Inf, ymin = Q1, ymax = Q3, fill = layers)) +
    geom_jitter(data = DAT, aes(x=.data[[feat]], y={{spatial_dist}}, col=layers), 
                width = 0.1, alpha = 0.7, size=.3) + #
    scale_fill_manual(values = col) + 
    scale_colour_manual(values = col) +
    # geom_smooth(data = filter(DAT, .data[[feat]] != 0), n=1000, aes(y={{spatial_dist}}, x=.data[[feat]], col=orig.ident)) + 
    guides(fill = guide_legend(override.aes = list(size=2), keyheight = .7, keywidth = .7)) +
    scale_y_reverse(expand = c(0, 0)) +
    #scale_x_continuous(expand = c(0, 0)) +
    {if(!(is.null(x_max))){xlim(-.5, x_max)}} +
    facets +
    my_theme +
    theme(plot.margin = unit(c(.2,0,0,.2), "lines"),
          #legend.box.margin=margin(0,0,0,0),
          legend.margin=margin(0,0,0,-5),
          panel.border = element_blank(),
          axis.line = element_line(),
          panel.grid.major = element_line(linewidth = 0.2),
          panel.grid.minor = element_line(linewidth = 0.1))
  return(dot)
}

col <- c("#E41A1C","#FF7F00","#C77CFF","#984EA3")
dot_epi <- map(cond_epi_DEGs, ~layer_dotplot.fun(DATA, .x, sp_dist_epi))
col <- c("#984EA3","#00BFC4","#00A9FF","#377EB8","#CD9600","#e0e067","#7CAE00","#FFFF33")
dot_sub <- map(cond_epi_DEGs, ~layer_dotplot.fun(DATA, .x, sp_dist_SubMuc, morf="SubMuc", clus=sub_clus))

############################
# COND EXPRESION ON TISSUE #
############################
# col <- RColorBrewer::brewer.pal(9,"PuRd")
# col <-  c("grey95", RColorBrewer::brewer.pal(9,"Reds"))
# col <- c("grey100","grey95", "mistyrose", "red", "dark red", "#870808", "black")
# col <- RColorBrewer::brewer.pal(9,"Purples")
col <- c("#EFEDF5", "#DADAEB", "#BCBDDC", "#9E9AC8", "#807DBA", "#6A51A3", "#54278F", "#3F007D") # Purples

# cond_epi_DEGs <- c("SAMD9", "GPRC5A", "TGM3", "KRT19", "PKP1")
tissue_epi <- map(cond_epi_DEGs, 
        ~plot_st_feat.fun( DATA,
                           geneid = .x,
                           zoom = "zoom",
                           col = col,
                           alpha = .9,
                           ncol = 2, 
                           save_space = F,
                           annot_line = .1,
                           img_alpha = 0,
                           point_size = .75)) 

tissue_sub <- map(cond_SubMuc_DEGs, 
        ~plot_st_feat.fun( DATA,
                           geneid = .x,
                           zoom = "zoom",
                           col = col,
                           alpha = .9,
                           ncol = 2, 
                           annot_line = .1,
                           img_alpha = 0,
                           point_size = .75)) 

# dev.new(width=6.75, height=5, noRStudioGD = TRUE)
wrap_plots(dot_epi[[1]],tissue_epi[[1]], ncol = 2, widths = c(.7, 1)) + plot_layout(guides = 'collect', axis_titles = "collect" )

#############################
# CONDITION DOTPLOT AS PDF #
#############################
# dev.new(width=(1.05*10)+1, height=2.5, noRStudioGD = TRUE)  
dot_epi <- wrap_plots(dot_epi, ncol = 5 ) + plot_layout(guides = 'collect', axis_titles = "collect")
ggsave("./Figures/08/cond_tissue_top_DEGs_epi.pdf", dot_epi, width = (1.4*5)+1, height = 2.5*ceiling(length(cond_epi_DEGs)/5), limitsize = FALSE )

dot_sub <- wrap_plots(dot_sub, ncol = 5 ) + plot_layout(guides = 'collect', axis_titles = "collect")
ggsave("./Figures/08/cond_tissue_top_DEGs_sub.pdf", dot_sub, width = (1.4*5)+1, height = 2.5*ceiling(length(cond_SubMuc_DEGs)/5), limitsize = FALSE )

# pdf("./Figures/08/condition_dotplot_DEGs.pdf",  width = (1.05*10)+1, height = 2.5*ceiling(length(feat_cond)/5))
# wrap_plots(dot, ncol = 5, nrow = ceiling(length(feat_cond)/5)) + plot_layout(guides = 'collect')
# dev.off()

#################################
# CONDITION TISSUE PLOTS AS PDF #
#################################
# dev.new(width=4, height=6, noRStudioGD = TRUE) 
# p[[1]]
ggsave("./Figures/08/condition_tissue_DEGs_epi.pdf", gridExtra::marrangeGrob(grobs = tissue_epi, nrow=1, ncol=1, top=NULL), width = 9, height = 3.5)
ggsave("./Figures/08/condition_tissue_DEGs_SubMuc.pdf", gridExtra::marrangeGrob(grobs = tissue_sub, nrow=1, ncol=1, top=NULL), width = 9, height = 3.5)

Plot cluster diff. gene expression as dotplot and on tissue

# dev.new(width=2, height=2.5, noRStudioGD = TRUE) 
# dev.new(width=(1.2*5)+1, height=2.5*ceiling(length(dot_fig)/5), noRStudioGD = TRUE) 
##############################
# CLUS DOTPLOT PER EPI LAYER #
##############################
col <- c("#E41A1C","#FF7F00","#C77CFF","#984EA3")
dot_epi <- map(cond_epi_DEGs, ~layer_dotplot.fun(DATA, .x, sp_dist_epi, facet = F))
col <- c("#984EA3","#00BFC4","#00A9FF","#377EB8","#CD9600","#e0e067","#7CAE00","#FFFF33")
dot_sub <- map(cond_epi_DEGs, ~layer_dotplot.fun(DATA, .x, sp_dist_SubMuc, morf="SubMuc", clus=sub_clus, facet = F))

############################
# CLUS EXPRESION ON TISSUE #
############################
tissue_epi <- map(clus_epi_DEGs[1], 
        ~plot_st_feat.fun( DATA,
                           geneid = .x,
                           zoom = "zoom",
                           col = col_feat,
                           alpha = .9,
                           ncol = 4, 
                           annot_line = .1,
                           img_alpha = 0,
                           point_size = .7))

tissue_sub <- map(clus_SubMuc_DEGs[1], 
        ~plot_st_feat.fun( DATA,
                           geneid = .x,
                           zoom = "zoom",
                           col = col_feat,
                           alpha = .9,
                           ncol = 4, 
                           annot_line = .1,
                           img_alpha = 0,
                           point_size = .7))

# dev.new(width=6.75, height=5, noRStudioGD = TRUE)
wrap_plots(dot_epi[[1]],tissue_epi[[1]], ncol = 2, widths = c(.5, 1)) + plot_layout(guides = 'collect', axis_titles = "collect" )

###########################
# CLUSTERS DOTPLOT AS PDF #
############################
dot_epi <- wrap_plots(dot_epi, ncol = 5) + plot_layout(guides = 'collect', axis_titles = "collect")
dot_sub <- wrap_plots(dot_sub, ncol = 5) + plot_layout(guides = 'collect', axis_titles = "collect")
ggsave("./Figures/08/clus_tissue_top_DEGs_epi.pdf", dot_epi, width = (1.2*5)+1, height = 2.5*ceiling(length(dot_epi)/5), limitsize = FALSE )
ggsave("./Figures/08/cond_tissue_top_DEGs_sub.pdf", dot_sub, width = (1.2*5)+1, height = 2.5*ceiling(length(cond_SubMuc_DEGs)/5), limitsize = FALSE )

# pdf("./Figures/08/cluster_dotplot_top_DEGs.pdf",  width = (1.2*5)+1, height = 2.5*ceiling(length(feat_clus)/5))
# wrap_plots(dot_fig, ncol = 5, nrow = ceiling(length(feat_clus)/5)) + plot_layout(guides = 'collect')
# dev.off()
#################################
# CLUSTERS TISSUE PLOTS AS PDF #
#################################
ggsave("./Figures/08/clus_tissue_top_DEGs.pdf", gridExtra::marrangeGrob(grobs = tissue_epi, nrow=1, ncol=1, top=NULL), width = 9, height = 3.5)
ggsave("./Figures/08/clus_tissue_top_DEGs.pdf", gridExtra::marrangeGrob(grobs = tissue_sub, nrow=1, ncol=1, top=NULL), width = 9, height = 3.5)